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This  instrumentation  grant  has  enabled  a  dramatic  enhancement  in  in-house  supercomputing 
power  at  ACMS.  The  SGI  0NYX2  22-processor  system  with  7  Gbytes  of  memory  was  replaced 
by  the  new  SGI  Altix  32-Itanium  processor  system  with  64  Gbytes  of  memory.  The  combination 
of  higher  speed  per  processor  (factor  of  6),  higher  interprocessor  bandwidth  and  much  larger 
accessible  shared  memory  (64  Gbytes)  has  meant  that  problems  that  were  previously 
inaccessible  are  now  well  within  reach.  Moreover,  this  system  offers  an  exciting  upgrade  path, 
initially  to  dual-core  and  eventually  to  multi-core  (16)  processors.  This  new  capability,  installed 
in  July  2004,  is  already  having  a  significant  impact  on  the  three  AFOSR  funded  projects: 

1 .  F496200210380  Novel  Designs  and  Coupling  Schemes  for  Affordable  High  Energy  Laser 
Modules 

2.  F4962003 10194  Computational  Nonlinear  Optics:  Femto  Atmospheric  Light  String 
Applications 

3.  FA95500410213  Fundamental  Modeling  and  Design  Strategies  in  Computational 
Photonics:  Applications  to  Lasercomm  through  clouds  and  Electro-optical/Nanophotonics 

In  addition,  the  added  leverage  provided  by  the  purchase  of  the  Altix  system  and  a  very  generous 
trade-in  value  on  the  old  ONYX2  system,  enabled  us  to  purchase  a  21 -node  (42 -CPU)  AMD 
Opteron  cluster  with  84  Gbytes  of  memory.  This  cluster,  which  has  a  high  speed  Infiniband 
interconnection  fabric,  took  some  time  to  stabilize  as  it  represented  the  first  such  cluster  system 
integration  attempt  involving  Silicon  Graphics  and  Atipa  engineers.  The  goal  was  to  seamlessly 
integrate  the  Altix  shared  memory  system  with  this  mixed  shared/distributed  memory 
environment  thereby  offering  unparalleled  supercomputing  performance.  The  full  system 
integration  involved  three  vendors:  Silicon  Graphics  (SGI  Altix),  Atipa  (Opteron  Cluster)  and 
Voltaire  (infiniband  interconnection  fabric  between  Cluster  nodes,  fast  local  node  Virtual 
Parallel  1  TB  disk  storage,  3  TB  disk  storage  and  the  connection  to  the  Altix).  The  fully 
integrated  mixed  supercomputing  environment  became  available  in  July  2005  with  the 
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installation  of  a  ROCKS  Redhat  Linux  queuing  system.  This  integrated  supercomputing 
laboratory  is  illustrated  in  the  figure  attached.  We  have  just  replaced  8  dual  Opteron  nodes  on  the 
cluster  with  8  dual  core,  dual  CPU  modules.  This  has  increased  the  number  of  available  cluster 
CPUs  further  by  16  processors  and  moreover,  provides  8  quad  processor  nodes  for  smaller  scale 
shared  memory  computing. 

F496200210380  Novel  Designs  and  Coupling  Schemes  for  Affordable  High  Energy 
Laser  Modules 

This  MRI  project  has  benefited  greatly  from  the  new  supercomputing  environment.  Large  scale 
simulations  of  double  clad  and  micro-structured  active  fibers  is  now  possible.  Record  powers 
with  single  frequency  operation  at  eye-safe  1.5  pm  wavelength  from  centimeter  long  Er/Yb- 
doped  phosphate  glass  doped  fibers.  New  large  core  micro-structured  phosphate  fiber  designs  are 
currently  being  investigated  as  a  pathway  to  further  power  scaling  by  shortening  further  the 
active  fiber  to  achieve  even  greater  longitudinal  mode  discrimination.  The  second  project  on 
power  scaling  of  semiconductor  VECSEL  structures  has  received  a  major  boost  from  the 
enhanced  supercomputing  performance  afforded  by  the  SGI  and  Cluster  systems.  Our  unique 
closed-loop  approach  to  designing  such  high  power,  high  brightness  laser  systems  requires 
supercomputing  at  many  different  levels.  First,  we  need  to  optimize  the  semiconductor  epi 
structure  for  operating  laser  active  layer  temperatures.  This  requires  microscopically  computing 
the  semiconductor  quantum  well  (QW)  optical  gain  and  refractive  index  spectra  over  the  relevant 
physically  accessible  landscape  of  carrier  density  and  temperature.  These  properties  together 
with  wafer-level  photoluminescence  spectra  are  proving  invaluable  in  providing  direct  feedback 
to  the  semiconductor  materials  growers.  Large  scale  simulation  including  optical  and  thermal 
transport  in  the  semiconductor  active  mirror  (sub-cavity),  free-space  propagation  to  the  external 
passive  mirror  and  full  thermal  transport  though  all  heat-spreading  and  heat-sinking  elements 
provides  a  major  computational  challenge.  We  currently  have  2  versions  of  such  a  simulator 
working  on  our  SGI  shared  memory  machine:  1)  A  full  scale  simulation  with  coarse-grained 
time  stepping  -this  provides  a  powerful  tool  for  designing  and  optimizing  the  whole  device  and, 
2)  Treatment  of  the  semiconductor  active  mirror+cooling  elements  as  a  lumped  active  mirror 
(with  inputs  provided  from  1))  and  a  full  spatio-temporal  simulation  of  the  transverse  and 
longitudinal  mode  structure.  The  latter  will  allow  us  to  optimize  the  active  VECSEL  for 
spectrally  narrow,  high  brightness  performance. 

F496200310194  Computational  Nonlinear  Optics:  Femto  Atmospheric  Light  String 
Applications 

The  main  effort  here  has  been  the  development  of  ultra-fast,  fully  vectorial,  vector  Maxwell 
solvers  for  applications  to  the  areas  of  atmospheric  light  string  propagation,  self-trapping  of 
pulsed  lasers  in  condensed  media,  white-light  supercontinuum  generationin  bulk  media  and 
micro-structured  fibers  and  nonlinear  X-waves.  These  codes  have  been  implemented  in  parallel 
on  the  SGI  Altix  machine  and  show  excellent  scaling  to  multiple  CPUs.  A  significant  number  of 
papers  have  been  published  based  on  these  solvers  and  utilizing  the  SGI  compute  engine.  These 
have  been  documented  in  annual  reports  for  this  grant.  A  major  new  thrust  of  the  research  is  to 
design  more  rigorous  physical  models  for  ultra-fast  light-matter  interaction  and  couple  these  to 
our  unidirectional  Maxwell  and  full  Maxwell  solvers.  All  computational  aspects  of  these  models 


are  extremely  challenging  and  we  are  collaborating  with  different  groups  worldwide  in  search  of 
computationally  efficient  but  physically  accurate  computational  models  and  schemes.  We 
envisage  taking  advantage  of  our  unique  mixed  shared/distributed  memory  computer  architecture 
in  tackling  some  of  these  emerging  problems.  All  efforts  here  are  particularly  relevant  to  the 
emerging  research  field  of  extreme  nonlinear  optics.  A  manuscript,  detailing  the  main  theoretical 
developments  of  this  work  is  attached. 

FA95500410213  Fundamental  Modeling  and  Design  Strategies  in  Computational 
Photonics:  Applications  to  Lasercomm  through  clouds  and  Electro-optical/Nanophotonics 

The  main  thrust  of  the  initial  phase  of  this  effort  has  been  in  researching  current  state-of-the-art 
theoretical  and  computational  methodologies  with  a  view  to  exploiting  our  progress  made  under 
F496200210380  on  the  development  of  high  power,  high  brightness,  spectrally  narrow  VECSEL 
sources.  Semiconductor  wafers  grown  under  the  current  project  are  designed  from  scratch  by  the 
ACMS  research  team.  Currently  we  are  exploring  whether  it  is  feasible  to  design  partially- 
coherent  lasercomm  sources  based  on  VECSEL  arrays. 

The  computational  nanophotonics  effort  relies  heavily  on  the  distributed  memory  cluster  for  3D 
vector  Maxwell  simulation  using  our  implementation  of  the  FDTD  time  domain  approach.  Our 
working  non-uniform  mesh  FDTD  Maxwell  solver  is  being  employed  to  investigate  light 
coupling  from  external  sources  into  photonic  crystal  and  nano-feature  structures.  We  are  also 
investigating  the  full  3D  confined  defect  modes  in  3D  photonic  crystals  with  a  view  to  isolating 
the  localized  evanescent  field  from  the  propagating  component.  The  role  of  surface  plasmon 
modes  in  enhancing  transmission  through  nanoscale  holes  in  metallic  films  is  also  being  studied. 
The  dramatically  increased  throughput  offered  by  our  supercomputing  facility  has  significantly 
enhanced  our  capability  to  solve  new  classes  of  computationally  intensive  problems. 


ACMS  Photonics  Supercomputing  Laboratory 


2  TB  PVFS  Scratch  Directory 
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Full  Vectorial,  Intense  Ultrashort  Pulse 
Propagators:  Derivation  and  Applications 


J.  V,  Moloney  and  M.  Kolesik 

Arizona  Center  for  Mathematical  Sciences  and  Optical  Sciences  Center, 

University  of  Arizona,  Tucson  AZ  85721 
jmlSacms . arizona . edu 

1.1  Introduction 

Rapid  progress  in  recent  years  in  the  development  of  high  power  ultrashort 
pulse  laser  systems  has  opened  up  a  whole  new  vista  of  applications  and  com¬ 
putational  challenges.  Amongst  those  applications  that  are  most  challenging 
from  a  computational  point  of  view  are  those  involving  explosive  critical  self- 
focusing  with  concomitant  explosive  growth  in  the  generated  light  spectrum. 
Moreover,  new  experimental  developments  in  the  field  of  extreme  nonlinear 
optics  will  require  more  rigorous  propagation  models  beyond  those  existing 
in  the  current  literature.  Specific  aplications  areas  chosen  for  illustration  in 
this  chapter  include  atmospheric  light  string  propagatio  and,  nonlinear  self¬ 
trapping  in  condensed  media.  These  examples  exhibit  rather  different  aspects 
of  intense  femtosecond  pulse  propagation  and  demonstrate  the  robustness 
and  flexibility  of  the  unidirectional  Maxwell  propagator.  A  novel  aspect  of 
our  approach  is  that  the  pulse  propagator  is  designed  to  faithfully  capture 
the  light-material  interaction  over  the  broad  spectral  landscape  of  relevance 
to  the  interaction.Moreover  the  model  provides  a  seamless  and  physically  self- 
consistent  means  of  deriving  the  many  ultrashort  pulse  propagation  equations 
presented  in  the  literature. 


1.2  Derivation  of  Unidirectional  Pulse  Propagation 
Equations 

1.2.1  Starting  point:  Maxwell  Equations 

In  this  Section  we  expand  the  discussion  presented  in  [1],  We  outline  the 
key  steps  in  deriving  a  physically  self-consistent  and  robust  ultrashort  pulse 
propagator  that  resolves  the  underlying  optical  carrier  wave  while  enabling 
propagation  over  many  meter  propagation  lengths.  Our  goal  is  to  retain  the 
full  rigor  of  Maxwell’s  equations  while  reducing  the  problem  complexity  by 
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constraining  the  model  to  unidirectional  propagation.  As  our  immediate  inter¬ 
est  is  in  very  short  intense  pulse  propagation  with  potentially  large  induced 
nonlinear  polarization,  we  will  need  to  accurately  capture  the  very  broad 
spectral  landscape  that  the  pulse  experiences  during  its  interaction  with  a 
host  dielectric  material.  In  many  cases,  spectral  superbroadening  is  such  that 
the  generated  bandwidth  far  exceeds  in  magnitude  the  underlying  carrier  fre¬ 
quency  i.e  Auj/ojo  »  1.  In  this  limit,  we  expect  the  Nonlinear  Schrpdinger 
Equation  (NLSE)  to  fail.  Many  attempts  have  been  made  to  derive  nonlin¬ 
ear  envelope  models  that  go  beyond  NLSE  and  we  will  discuss  most  of  these 
below  when  we  show  explicitly  how  each  can  be  seamlessly  derived  from  our 
unidirectional  pulse  propagation  equation  (UPPE). 


Time-propagated  and  space-propagated  equations 

Most  of  the  pulse  propagation  problems  in  nonlinear  optics  are  solved  in  one 
of  two  formulations:  Either  one  has  an  initial  condition  (electric  and  magnetic 
fields)  specified  in  all  space,  and  the  evolution  is  calculated  along  the  time 
axis,  or  the  initial  condition  is  given  as  a  function  of  local  pulse  time  and 
transverse  (w.r.t.  propagation  direction)  coordinates,  and  the  numerical  evo¬ 
lution  proceeds  along  the  propagation  axis.  We  refer  to  these  cases  as  time- 
and  2-propagated  equations. 

The  z-propagated  approach  is  much  more  common  in  nonlinear  optics 
simulations  based  on  envelope  equations,  often  related  to  NLS.  The  time- 
propagated  approach  is  on  the  other  hand  common  for  solvers  based  on  direct 
integration  of  Maxwell’s  equations. 

Due  to  space  limitations,  we  focus  in  this  chapter  on  the  z-UPPE.  As 
discussed  in  [1]  in  more  detail,  the  time-propagated  versions  of  UPPE  are 
more  suitable  for  tight-focusing  scenarios  when  non-paraxial  effects  start  to 
play  a  role.  The  2-propagated  equations  are  easier  to  use  in  situations  that 
allow  to  neglect  the  longitudinal  field  components  as  sources  of  nonlinear 
material  responses. 

The  2-propagated  approach  is  more  common  in  nonlinear  optics  and  we 
start  by  deriving  it  explicitly,  then  discuss  briefly  its  numerical  implementa¬ 
tion,  relation  to  other  envelope  models  and  provide  some  illustrative  exam¬ 
ples  of  its  implementation.  We  refer  the  reader  to  [1]  for  details  concerning 
the  details  of  numerical  treatment  of  the  nonlinear  material  response  in  both 
t-UPPE  and  z-UPPE  cases. 


1.2.2  Starting  point:  Maxwell  Equations 

As  a  first  step  in  derivation  of  various  versions  of  UPPE,  we  derive  an  exact 
coupled-modes  system  of  equations.  This  is  a  well-known  textbook  formula 
that  can  be  found  in  the  literature  in  several  different  forms.  To  keep  our 
derivation  self-contained,  we  re-derive  the  starting  formula  using  a  standard 
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approach  based  on  orthogonality  relations  for  modes  of  an  electromagnetic 
field. 

Electromagnetic  fields  of  a  light  pulse  propagating  along  the  z  axis  can  be 
expanded  into  modal  contributions  that  reflect  the  geometry  of  the  waveguide 
(we  consider  homogeneous  medium  a  special  case  of  the  latter). 

E(x,y,z,t)  =  ^Arn{u,z)£rn{u,x,y)ei0m{“)z~lU3t 

m,(jj 

H(x,y,z,t)  =  (1.1) 

m,L> 


Here,  m  labels  all  transverse  modes,  and  an  initial  condition  Am(uj,  z  =  0)  is 
supposed  to  be  given  or  calculated  from  the  field  values  at  2  =  0.  Note  that 
the  above  expansion  is  valid  for  the  transverse  components  only. 

To  save  space,  the  following  short-hand  notation  will  be  used  below 

£m  =  Cm(«,*,y)e<ft-<w>-fcrt 

Hm  3  .  (1.2) 

The  starting  point  is  Maxwell’s  equations.  We  consider  a  non-magnetic 
medium  (p  =  po)  with  a  linear  permittivity  e(u>,  x,  y)  that  doesn’t  depend  on 
the  propagation  coordinate  2. 

j  +  dtP  +  eodte  *  E  =  V  x  H 

-li0dtH  =  Vx£  (1.3) 

where  the  star  represents  a  convolution. 

First,  we  scalar-multiply  the  Maxwell’s  equations  by  the  complex  conjugate 
modal  fields 


+  dtP)  +€0 £*m-dte  *  E  =  £*m.V  x  H 

-y0H*m.dtH  =  H*m.V  x  E  .  (1.4) 

Using  the  formula  6.(V  x  a)  =  V.(a  x  6)  +  a.(V  x  b),  we  transform  both 
RHS  to  obtain 

£*m.(j  +  dtP)  +  e0£*m.dte  *  E  =  V.[H  x  £*m]  +  H.[V  x  £*m) 

-HoH*m.dtH  =  V.[£?  X  +  E.[V  X  H*m}  .  (1.5) 

Since  modal  fields  themselves  satisfy  the  Maxwell’s  equations 

V  x  £*m  —  -po dtU'm 

VxH*m  =  e  0dte*£*m,  (1.6) 

the  above  equations  can  be  written  as 

£*m-U  +  dtP)  +  eo£*m.dtt  *  E  =  V.[H  x  £*m]  - 

-po H*m.dtH  =  V.[E  x  H*m]  +  e0E.dte  *  £*m 


(1.7) 
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We  subtract  the  two  equations  and  collect  terms  that  constitute  full  time 
derivatives 

£*m-U  +  dtP )  +  dt[e0S*m.e *  E )  =  V.fff  x  £*J  -  dt\»0 -V.[Ex  H*m]  . 

(1.8) 

Now  integrate  over  the  whole  xyt  domain.  Note  that  all  terms  except  the 
first  and  the  dz  (implicit  in  V.)  are  derivatives  that  give  rise  to  “surface 
terms”  after  integration  over  x,  y,  t.  Since  we  consider  spatially  and  temporally 
localized  pulse-like  solutions,  these  terms  vanish.  The  only  surviving  derivative 
will  be  dz  from  V.: 

Js-m.(j  +  d,P)ixdyit  =  a.Jz.lH  x£*m]dxdydt-dz  J  z.[E  x  7i.*m]dxdydt 

(1.9) 

Here  and  in  what  follows,  t  integrations  are  understood  like  this:  f  dt  = 
j,  J^t/2  ^  where  T  is  a  large  normalization  “volume,”  and  integrals  over  x,  y 
are  understood  in  a  similar  way.  Because  only  transverse  field  components 
enter  above  equation,  we  can  use  our  modal  expansion  here: 

f£*m.(j  +  dtP)dxdydt  = 

&/*•[£„, fl^(^,z)?in(/2)  x  £Uu))eiMn)z~inte~i0m(w)z+iutdxdydt 
-dz  f  z.[£n' «  An(fi,z)£n(n)  x  WMje^^-^e-^^+^dxdydt  . 

(1.10) 

Integration  over  time  gives  a  Kronecker  delta  between  angular  frequencies, 
Snu, ,  which  in  turn  reduces  the  sum  over  i?: 

!  £*m-{j  +  dtP)dxdydt  = 

dz  J  2-[£n  An{w,  z)Hn(u,  x ,  y )  x  £*m{ui,  x ,  j/)]ei/3”(“’^e_i/J”(a,)zdxd?/ 

-dz  f  z.[EnAn(u,z)£n(u,x,y)  x  x,  dxdy  . 


Collecting  like  terms  results  in  an  equation 


[  £m-U  +  9tP)dxdj/dt  =  dz  ^2  An(u,  z)eil3n(^ze 

n 

I**1  n(w,  x,  y)  x  £*m{u>,  x,  y)  -  £n{u,  x,  y)  x  W^(w,  x,  2/)]dxd?/  .  (1.12) 


where  the  general  orthogonality  relation 


J  z.[£m  x  "H.n  T~Cm  x  £n\  dxdy  —  2 Sm>nNm(u})  (1.13) 


can  be  used  to  reduce  the  sum  over  modes 


J +  dtP)dxdydt  =  -dzy}TAn{u>1z)ei0''^ze  if3m{-u'>z28rn%nNm{u)  , 
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and  finally  obtain  the  evolution  equation  for  the  expansion  coefficients: 
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^  <*+T/2  /•+Y/2  p+X/2 

-xytJ_t/1  *J_r/2  *vj_xn 

e-iM")M£*mi^Xjy)'{ j(x,y,t)  +  dtP(x,y,t)}  (1.15) 
This  is  the  starting  point  for  the  z-propagated  unidirectional  equations. 

1.2.3  Z-propagator  UPPE  for  homogeneous  medium 

In  the  following,  Eq.  (1.15)  is  specialized  for  the  case  of  a  homogeneous 
medium.  Field  modes  are  plane  waves  labeled  by  transverse  wavenumbers 
kx,  ky,  a  polarization  index  s  =  1, 2,  and  a  ±  sign  depending  on  the  direction 
of  propagation  along  z: 

m  =  hx,ky,s,±  .  (1-16) 

The  frequency  and  wavenumber  dependent  propagation  “constant”  is 

&x,*»,«,±(w)  =  kz(u>,  kx ,  ky)  =  ^Jui2e(uj)/c2  -k2  -k2  ,  (1.17) 

and  the  modal  field  amplitudes  are 

£kx,ky,s,±  =  es  exp  [ ikxx  +  ikyy  ±  ikz(u>,  kx,  ky )]  (1.18) 

'kLkx,ky,s,±  —  k  X  £kx,ky,w,s,±  ■  (1-19) 

flQUJ 

Here,  ea=i,2  are  unit  polarization  vectors  normal  to  the  wave-vector 


k  =  {fc*,  ky,  kz  =  y'w 2e(w)/c 2  -  k2  -  k2}  .  (1.20) 

From  these  formulas,  it  is  straightforward  to  calculate  the  modal  normal¬ 
ization  constant 

2 Nkx}kv,St±(w)  =  J  z.[£m  xTi^-HmX  £*J  dxdy  = 

1  1 

2z.[es  x  (k  x  es)] - =  ±2 kz(u>,  kx,  kz) -  (1-21) 


NkXlky,S,±(v)  =  ± 


kz  (tu,  kx,  kz) 


Now  we  can  insert  expressions  for  the  modal  fields  and  normalization  con¬ 
stant  into  coupled  mode  equation  Eq.  (1.15) 

dzAk  k  s  +(w,z)  =  eTikzZ  f  dxdydtei(^-kxx-kvy) 

z  k*'ky's’+ V  '*>  2 kz  J  LxLyT 

es-[j{x,  y,  z,  t )  +  dtP{x,  y,  z,  t)]  (1.23) 
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The  above  integral  is  nothing  but  spatial  and  temporal  Fourier  transform,  so 
one  can  write  it  in  the  spectral  domain  as 


dzAk, 


,ky 


_(w,z)  = 


U) 


2eo  (?k- 


■e  lk*zea.[iuPkx,ky(u,z)  -  jkx,ky(u,z)]  (.1.24) 


This  is  the  propagation  equation  that  is  actually  solved  numerically,  but  let  us 
express  it  in  terms  of  field  rather  than  in  terms  of  modal  expansion  coefficients. 
From  a  modal  expansion,  the  transverse  part  of  the  electric  field  is 


El 


=  E  ejAkx,kvtSt+^z)eik^ky^  ,  (1.25) 

8=1,2 


and  its  z  derivative  reads 


dzEl,ky,+{u,z)  =  ikz(kx,kyt<j)El'kv'+{u,z)+ 

E  ejdzAkxiky,s,+(u,z)eik>tk*>k»^z  (1.26) 

8=1,2 


Using  Eq.  (1.24)  we  obtain  the  homogeneous  medium  full-vectorial  UPPE 

d*El,kv,+{w'z)  =  ikzEkx<kv,+(u},z)+ 

•  2 

£ -[^k;Pk*’k>'z)  ~  2^kJk-k>'z)]  (L27) 

5=1,2 

This  is  an  exact  system  of  equations  that  describes  the  evolution  of  modal  am¬ 
plitudes  along  the  z-axis  for  the  forward  propagating  field.  A  similar  equation 
holds  for  the  backward  propagating  component,  of  course. 

Because  the  nonlinear  polarization  in  this  equation  results  as  a  response 
to  the  complete  field,  this  equation  can’t  be  used  to  calculate  the  forward 
field  only.  The  equation  becomes  “unidirectional”  only  after  resorting  to  the 
following  approximation 


P(E),j(E)->P(Ef),j(Ef)  (1.28) 

In  other  words,  to  obtain  a  closed  system  to  solve  numerically,  we  must  re¬ 
quire  that  the  nonlinear  polarization  is  well  approximated  by  the  nonlinear 
polarization  calculated  only  from  the  forward  propagating  field.  This  means 
that  the  equation  is  only  applicable  when  the  back-reflected  portion  of  the 
field  is  so  small  that  its  contribution  to  the  nonlinearity  can  be  neglected. 


1.2.4  Z-propagated  UPPE,  simplified,  most  practiced  version 

Eq.  1.27,  with  nonlinear  polarization  approximated  by  Eq.  1.28  can  easily  be¬ 
come  a  rather  large  system  to  solve  numerically  for  typical  problems  encoun¬ 
tered  in  the  femtosecond  pulse  propagation  area.  Fortunately,  in  most  cases 
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the  transverse  dimensions  of  the  resulting  field  remain  relatively  large  com¬ 
pare  to  the  wavelength,  and  further  approximations  are  possible.  Concretely, 
one  can  neglect  the  z  components  of  the  field  and  polarization  vectors.  In  such 
a  situation  the  sum  over  polarization  vectors  reduces  approximately  to  unity 

£  es'es  «  1  .  (1-29) 

a=l,2 

and  the  full  UPPE  simplifies  into  an  equation  for  the  transverse  component  (s) 
of  the  field 

2Cj2  (jj 

dzEkx,ky(u,  z)  =  ikzEkx:ky(u,z)  +  2£oC2jfc  (w’ z)  ~  2g0 c2k  jkx’k«(u>,z')  ’ 
kz  =  1v/w2e(w)/c2  -  A:2  -  fc2  .  (1.30) 

This  is  the  most  useful  form  for  practical  calculation,  and  is  therefore  called 
simply  UPPE  in  the  following. 


1.2.5  Nonlinear  material  response 

In  most  cases,  the  propagation  equations  discussed  in  this  chapter  do  not 
require  a  specific  form  of  material  response.  However,  for  the  sake  of  concrete¬ 
ness,  as  well  as  for  discussion  of  numerical  methods,  we  want  to  describe  a 
generic  model  of  nonlinear  material  response.  We  consider  a  nonmagnetic,  dis¬ 
persive  medium  with  relative  permittivity  e  that  is  a  function  of  the  transverse 
coordinates  x,  y  and  of  the  angular  frequency  ui 

e  =  e(oj,x,y)  ,  n  =  no  ■  (1.31) 

This  medium  specification  includes  any  dispersive  homogeneous  medium  such 
as  air  or  water  as  well  as  structured  fiber-like  media  such  as  photonic,  mi- 
crostructured  and  tapered  optical  fibers. 

Nonlinear  effects  are  usually  described  in  terms  of  polarization  P  through 
the  material  constitutive  relation: 

D  =  e0e  *  E  +  P  .  (1.32) 

The  star  in  this  formula  represents  a  convolution  integral  with  e  being  the 
linear  response  function  corresponding  to  the  frequency  dependent  e(uj,x,y). 
The  non-linear  polarization  is  an  “arbitrary”  function  of  the  electric  field 
P  =  P(E).  We  will  also  include  a  current  density  that  is  driven  by  the 
optical  field 

3  =  3(E)  (1.33) 

to  describe  interactions  with  plasma  generated  by  the  high-intensity  optical 
pulse. 
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The  main  physical  effects  that  influence  propagation  of  ultrashort,  high- 
power  light  pulses  in  nonlinear  dispersive  media  include  the  optical  Kerr  and 
stimulated  Raman  effects,  free-electron  generation,  defocusing  by  the  gen¬ 
erated  plasma  and  losses  caused  by  avalanche  and  multiphoton  ionization 
(MPI).  With  minor  modifications,  models  including  these  effects  can  be  used 
for  description  of  ultra-short  optical  pulses  propagation  in  gases  [2-18],  con¬ 
densed  bulk  media  [19-23]  and  in  conventional,  microstructured,  and  tapered 
fibers  [24-26]  as  well  as  in  ultra-thin  silica  “wires”  [27]. 

The  optical  Kerr  and  stimulated  Raman  effects  cause  a  local  modification 
of  the  optical  susceptibility 

P  =  e0AXE  (1.34) 

that  responds  to  the  history  of  the  light  intensity  I : 


poo 

Ax  =  2 nbn2  (1  -/)/  +  //  H(T)I(t  -  t)<1t 

Jo 


(1.35) 


Here,  /  is  the  fraction  of  the  delayed  nonlinear  response,  and  7 Z  is  the  mem¬ 
ory  function  of  the  stimulated  Raman  effect.  Parameterization  by  TZ(t)  ~ 
sin(J7r)e-rT  is  often  sufficient  for  ultrashort  pulses  [28].  This  simple  formula 
has  the  advantage  of  easy  implementation  that  avoids  explicit  calculation  of 
the  convolution  integral.  Often,  an  even  simpler,  exponential  memory  func¬ 
tion  is  used,  TL(t)  ~  e~rT  in  simulations  (see  e.g.  [29]).  If  the  real  memory 
function  is  sufficiently  complex,  a  numerical  convolution  approach  must  be 
used  to  calculate  the  convolution.  This  is  e.g.  the  case  in  silica  [30]. 

Of  course,  E.  1.34  neglects  the  dependence  of  the  Kerr  effect  on  wavelength. 
Although  Ax  may  exhibit  a  finite  memory  due  to  the  Raman  contribution, 
it  acts  on  the  instantaneous  value  of  E  only.  This  is  in  part  due  to  only 
rather  limited  data  available  on  frequency  dependence  of  the  nonlinear  coef¬ 
ficients  ri2  (see  Ref.  [31]  for  silica),  but  it  also  simplifies  practical  calculations 
considerably.  Consequenly,  the  “background”  index  of  refraction  n  b  can  be 
approximated  by  a  constant  value  taken  at  the  central  frequency  of  the  initial 
pulse. 

Because  of  the  potentially  high  intensities  accuring  in  femtosecond  pulses, 
free  electrons  are  generated  by  MPI  and  avalanche  mechanisms.  Then  it  is 
necessary  to  account  for  the  response  of  the  optical  field  to  the  presence  of  a 
dilute  plasma.  Since  the  relevant  times  scales  are  so  short,  plasma  diffusion  and 
ion  motion  are  neglected,  and  the  free-electron  density  p  is  usually  obtained 
as  a  solution  to  an  equation  of  the  following  form  [12, 13,28] 

dtp  =  alp  +  &(/)  —  cp2  .  (1.36) 


Here,  I  is  the  light  intensity,  a  parameterizes  the  avalanche  free-electron  gen¬ 
eration,  and  b(I)  represents  the  Multi  Photon  Ionization  (MPI)  rate  that  is 
a  highly  nonlinear  function  of  the  intensity.  The  last  term  describes  plasma 
recombination. 
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We  assume  that  the  collective  electron  velocity  v  responds  to  the  optical 
field  and  the  total  current  density  is  governed  by  the  following  simple  equation 
(see  e.g.  Ref.  [32]) 

im  =  ~ p(t)E(t)  -  j(t)/rc  (1.37) 

ut  ttIq 

where  rc  is  the  mean  time  between  collisions  experienced  by  electrons.  This 
equation  is  solved  together  with  (1.36)  to  capture  effects  of  the  plasma  on  the 
propagation  of  the  optical  field,  namely  defocusing  due  to  plasma  and  plasma 
induced  losses. 

Losses  caused  by  multiphoton  ionization  are  usually  incorporated  as  either 
an  equivalent  current  (see  e.g  [32, 33])  or  an  imaginary  susceptibility  contri¬ 
bution  that  extracts  from  the  field  the  energy  needed  for  the  free-electron 
generation. 


1.2.6  Numerical  algorithms 

The  general  structure  of  both  t-  and  z-UPPE  propagation  equations  is  quite 
similar.  A  UPPE  equation  actually  comprises  a  large  system  of  ordinary  dif¬ 
ferential  equation  with  a  free  propagation  part,  and  a  nonlinear  coupling. 

An  important  feature  that  affects  the  numerical  solution  strategy  is  that 
these  equations  are  written  in  the  spectral  space,  either  in  the  three  dimen¬ 
sional  space  of  wave-vectors  (f-propagated  UPPE)  or  in  a  two-dimensional 
space  of  transverse  wave-vectors  plus  a  one  dimensional  angular-frequency 
space  (z-propagated  UPPE).  At  the  same  time,  the  nonlinear  material  re¬ 
sponse  must  be  calculated  in  the  real-space  representation.  Consequently,  a 
good  implementation  of  Fast  Fourier  Transform  is  essential  for  a  UPPE  solver. 

Due  to  the  nature  of  FFT,  or  spectral  transforms  in  general,  parallelization 
of  the  solver  is  more  suited  to  a  shared  memory  architecture  than  to  the 
Message  Passing  Interface  approach.  Further,  since  the  most  time-consuming 
portion  of  the  code  deals  with  calculating  the  RHS  of  the  equation,  we  choose 
to  apply  a  single  threaded  library  for  ODE  solvers  (gsl  in  our  case)  an  only 
parallelize  the  calculation  of  derivatives  needed  in  the  ODE  solver. 

For  concreteness,  we  describe  solution  of  the  z-propagated  UPPE  equation: 

dzEkx,ky{u,  z)  =  ikzEkx,ky{w,z)  +  2eoC2/j  ~  2e0 c2k  ^kx'kv^u),z^  1 

kz  =  yj u)2 e{u>) / c2  -k2-k2  .  (1.38) 

Suppose  we  have  a  solution  at  z  =  0,  and  want  to  propagate  it  over  a  distance 
corresponding  to  an  integration  step  Az.  Although  we  wrote  the  equation  in 
terms  of  field,  it  is  more  natural  to  solve  numerically  for  the  modal  expansion 
coefficients  to  eliminate  fast  oscillatory  terms.  We  therefore  express  the  field 
through  the  expansion  coefficients  A  which  are  actually  the  “native  solver 
variables:” 
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Ekxiky(u>,z)  =  Akx,ky(u,z)eik^k*'k»>  (1.39) 

to  transform  the  equation  into 


dzAkx<ky(uj,z)  =  e-'k*z[+ 


2eo  c2k: 


■PkX,ky(.U,Z)- 


UJ 


2  enc2k: 


-jkx,ky(oJ,z)llAQ) 


Depending  on  the  concrete  ODE  scheme  one  chooses  to  use,  the  ODE  solver 
may  require  one  to  calculate  the  derivative  (RHS)  at  various  values  of  2.  Non¬ 
linear  polarization  and  current  are  of  course  functions  of  the  electric  field 
which  in  turn  is  calculated  for  any  given  z  from  (1.39).  That  is  nothing  but 
applying  a  linear  propagator  (in  the  spectral  representation)  to  a  given  field. 
Thus,  to  evaluate  the  RHS  in  (1.40)  given  Akxtky(u},  z),  one  first  applies  (1.39) 
to  obtain  electric  field  in  the  spectral-space  representation.  Subsequently,  a 
real-space  representation  is  obtained  by  Fourier  transform.  ( Note  that  this  is  a 
“global”  operation  which  forces  each  parallel  execution  thread  to  access  large 
amount  of  data  “owned”  by  all  other  threads.  This  is  the  main  reason  why  the 
UPPE  solvers  are  easier  to  implement  on  a  shared  memory  machine.)  Having 
the  electric  field  in  real  space,  one  calculates  the  material  response  using  the 
model  described  in  the  previous  Section.  This  is  done  by  integrating  the  ma¬ 
terial  equations  along  the  time  dimension  independently  for  each  point  in  the 
transverse  x,  y  plane.  Consequently,  this  portion  of  the  code  is  straightforward 
to  parallelize.  As  a  result,  one  obtains  the  nonlinear  polarization  and  current 
in  the  real-space  representation.  This  is  then  Fourier  transformed  back  into 
spectral  space  to  finally  evaluate  the  RHS  in  the  propagation  equation  (1.40). 
An  ODE  solver  calls  this  procedure  as  many  times  as  it  needs  to  perform  a 
single  integration  step  and  thus  calculates  Akx<ky(u,  z  +  Az). 

A  final  operation  that  finishes  one  solver  step  is  to  calculate  the  field  in 
the  frame  moving  with  a  an  appropriate  group  velocity  so  that  the  pulse  is 
kept  in  the  center  of  the  computational  domain.  This  is  achieved  by  applying 
a  modified  linear  propagator 


Akx,ky  (w,  z  +  Az)  -»  Akxiky (w,  2  +  Az)exp{i{kz{u ,  kx,  ky)  -  - — —}Az} 

V9 

(1.41) 

The  group  velocity  that  enters  above  is  usually  calculated  from  the  linear 
dispersion  relation 

—  =  dukz(u,kx  =  0,ky  =  0)|W=WR  (1-42) 

V9 

at  the  frequency  that  corresponds  to  the  maximum  of  the  pulse  spectral  power 
(pulse  center  frequency).  Let  us  emphasize  that  this  last  operation  in  just 
a  change  of  coordinates,  so  the  introduction  of  the  reference  frequency  at 
this  point  doesn’t  change  the  light  propagation  in  any  way.  However,  it  is 
important  to  include  the  moving  frame,  because  it  makes  it  easier  and  more 
efficient  to  apply  absorbing  boundary  layers  in  the  computational  domain. 
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Without  a  moving  window,  the  computational  domain  would  we  periodically 
wrapped  into  itself  in  the  time  direction.  In  order  to  apply  artificial  absorption 
at  the  boundary  and  also  for  integrating  the  material  equations,  the  latter 
would  need  to  be  identified  at  each  step. 

Thus,  since  the  UPPE  is  actually  a  big  ODE  system,  the  solver  imple¬ 
mentation  is  simple  in  principle,  though  it  involves  significant  “auxiliary” 
calculations  described  above  to  supply  the  ODE  solver  with  the  routine  to 
calculate  the  derivatives. 

A  final  remark  in  this  Section  concerns  axially  symmetric  problems.  We 
usually  treat  these  in  radial  coordinates  and  apply  a  numerical  Hankel  trans¬ 
form  instead  of  the  Fourier  transform.  This  is  a  “slow”  transform  with  a  dense 
matrix,  but  due  to  the  relatively  small  computational  domain  radially  sym¬ 
metric  problems  require,  this  is  not  a  big  problem.  Alternatively,  one  could 
treat  such  situations  by  finite  differencing  in  the  radial  dimension,  but  it  would 
mean  accepting  additional  (paraxial)  approximation,  and  would  introduce  ar¬ 
tificial  numerical  dispersion  into  the  algorithm. 


1.3  General  method  for  derivation  of  various 
propagation  equations  from  UPPE 

Several  types  of  unidirectional  propagation  equation  are  widely  used  in  the 
nonlinear  optics  literature.  The  most  important  examples  are  Non-Linear 
Schrodinger  (NLS)  equation  [34],  Nonlinear  Envelope  Equation  [35]  (NEE), 
the  First-Order  Propagation  equation  [33]  (FOP),  Forward  Maxwell’s  equa¬ 
tion  [36]  (FME),  and  several  other  equations  that  are  closely  related  to  these. 
The  derivations  found  in  the  literature  differ  from  equation  to  equation,  and 
in  some  cases  the  physical  meaning  of  the  required  approximations  may  not 
be  readily  evident  due  to  a  number  of  neglected  terms. 

In  this  section,  we  provide  a  unified  approach  that  will  subsequently  be 
used  to  derive  several  of  the  light-pulse  propagation  equations  found  in  the 
literature.  The  main  benefit  of  re-deriving  known  equations  from  a  common 
starting  point,  namely  UPPE,  using  the  same  method,  is  that  it  allows  us 
to  compare  the  physical  assumptions  and  approximation  underlying  different 
equations.  It  also  reveals  relations  between  the  equation  which  may  not  be  ob¬ 
vious  either  because  of  their  apparently  different  form,  or  because  of  different 
methods  used  in  the  original  derivations. 

It  is  instructive  to  break  the  derivation  procedure  into  several  steps.  As 
a  first  step,  we  adopt  a  scalar,  one  component  approximation  and  write  the 
UPPE  in  the  following  form: 

dzEk„,kv{u,z)  =iKEkx,kv{u>,z)  +  iQPks"kv(u,z)  (1.43) 


where 


K(kx,ky,w)  =  ^/w2e(u;)/c2  -  fc2  -  k\ 


(1.44) 
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is  the  linear  field  propagator  in  the  spectral  representation,  and 


J  ky  ,  CU) 


2 e0c2^u;2e(uj)/c2  -k%-k$ 


(1.45) 


will  be  called  nonlinear  coupling  term. 

In  the  second  step,  we  replace  K  and  Q  by  suitable  approximations.  In 
most  cases  they  are  nothing  but  Taylor  expansions  in  frequency  and  in  trans¬ 
verse  wavenumbers. 

To  obtain  envelope  equations,  one  expresses  the  field  in  terms  of  an  enve¬ 
lope  by  factoring  out  the  carrier  wave  at  a  chosen  reference  angular  frequency 
wr  with  the  corresponding  wave- vector  k r  =  ^(OjOjOJr): 

E(x,y,z,t)  =  A{x,y,z,t)e*k**-^  (1.46) 


A  similar  factorization  is  of  course  introduced  for  the  nonlinear  polarization 
P(x,y,z,t )  as  well. 

Step  three  consist  in  transforming  the  equation  from  the  spectral-  to  the 
real-space  representation.  Mathematically,  this  is  nothing  but  a  Fourier  trans¬ 
form  that  results  in  the  following  standard  rules  for  differential  operators: 

(lj  -  wr)  ->  idt  ikx  — >  dx  iky  — *  dy  dz  — ►  ik(ujR)  +  dz  (1-47) 

Finally,  in  most  cases  we  also  transform  to  a  frame  moving  with  a  suit¬ 
able  group  velocity  such  that  the  pulse  remains  close  to  the  center  of  the 
computational  domain. 


1.3.1  Derivation  of  Non-Linear  Schrodinger  Equation  from  UPPE 

The  Nonlinear  Schrodinger  Equation  (see  [34]  for  applications  in  optics)  (NLS) 
is  a  prototype  propagation  equation  in  nonlinear  optics.  This  is  also  the  sim¬ 
plest  case  suitable  to  illustrate  the  general  procedure  outlined  above.  One 
characteristic  feature  of  NLS  and  of  other  envelope  type  equations  is  the 
presence  of  a  reference  frequency.  Usually,  one  chooses  the  reference  angular 
frequency  wr  as  the  central  frequency  of  the  initial  pulse,  but  this  is  not  nec¬ 
essary.  Actually  it  is  useful  to  keep  in  mind  that  wr  is  to  a  certain  extent  a 
free  parameter,  and  that  the  obtained  results  must  be  almost  independent  of 
its  concrete  choice.  If  a  numerical  simulation  turns  out  to  be  sensitive  to  the 
choice  of  wr,  it  means  that  an  envelope  equation  is  being  used  outside  of  its 
region  of  validity. 

Following  the  general  procedure,  we  replace  the  K  and  Q  “coefficients” 
with  appropriate  approximations.  We  denote  by  A;r  =  k(u) r)  the  reference 
wavenumber  corresponding  to  the  chosen  reference  frequency  wr,  and  take 

K(kx,  ky,u>)  =  y/w2e(w)/c2  -k*-k% 
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kR  +  vg  x(uj  -  wR)  +  ^-(w  ~  WR)2  “  +  ^y)^1-48) 


This  is  a  second-order  Taylor  expansion  in  u)  —  wr  and  in  kx,ky. 

In  the  nonlinear  coupling  coefficient,  we  neglect  all  variable  dependencies 
and  take  its  value  at  the  reference  frequency  and  zero  transverse  wavenumbers: 


Qi^kxj  ky^Lj'j  - 


or 


wr 


2e0 Cyw2£(u)/c2  ~k2x-k2y  2e0 n(ojK)c 


(1.49) 


For  simplicity,  in  NLS  we  only  account  for  the  instantaneous  optical  Kerr 
effect,  and  write  the  nonlinear  polarization  envelope  as 


V  =  2eon(u)R)n2lA 


(1.50) 


Inserting  the  above  expressions  into  (1.43,1.46)  we  obtain 
dzA  =  +wJ1(uj-ur)A  +^(u-u}R)2A--^(kl+kl)A+t~n2IA  (1.51) 

It  is  customary  to  normalize  the  envelope  amplitude  such  that  \A\2  =  /.  Using 
rules  (1.47)  we  finally  obtain  the  NLS  equation: 

(dz  +  v~ldt)A  =  -  l\dttA  +  ^n2\A\2A  (1.52) 

The  above  derivation  shows  explicitly  what  approximations  need  to  be  adopted 
to  obtain  NLS:  Approximating  K  to  second  order  in  frequency  and  transverse 
wavenumber  amounts  to  the  paraxial,  and  quasi-monochromatic  approxima¬ 
tions  for  the  linear  wave  propagation.  The  approximation  in  the  nonlinear 
coupling  Q  also  requires  a  narrow  spectrum  in  order  to  be  able  to  represent 
Q  by  a  constant. 


1.3.2  Nonlinear  Envelope  Equation 

The  Nonlinear  Envelope  Equation  [35]  is  a  paraxial  equation  with  some  addi¬ 
tional  approximations  related  to  chromatic  dispersion.  This  equation  appears 
to  be  extremely  close  to  the  paraxial  version  of  UPPE. 

Once  again  we  follow  the  general  procedure  and  approximate  the  linear 
propagator  by  its  paraxial  version: 

K(kx,  ky,uj)  =  y/w2e(w)/c2  -  k%-k 2  «  +k(w)  -  2a,nf(^i'R)  +  kV>  t1’53) 

This  is  essentially  the  second-order  (paraxial)  Taylor  expansion  in  transverse 
wavenumbers  with  only  minor  additional  approximation.  Namely,  we  replaced 


14 


J.  V.  Moloney,  M.  Kolesik:  Ultrashort  Pulse  Propagators 


nb(oj)  — »  nb(v r)  in  the  denominator  of  the  diffraction  term,  and  thus  partly 
neglected  the  chromatic  dispersion. 

Further,  the  first  term  in  the  above  approximation,  which  is  an  exact 
propagation  constant  for  a  plane  wave  propagating  along  the  z  axis,  is  re¬ 
expressed  as  a  sum  of  its  two  lowest-order  Taylor  expansion  terms  plus  the 
rest: 

k(u>)  =  k(u>R )  -I-  v~ 1{u  -  wr)  -1-  D(u  -  wr)  (1-54) 


where 


D{u  -  wr) 


(oj  -  qjR)n 
n! 


(1.55) 


This  is  formally  exact  and  can  be  practically  implemented  in  the  spectral  do¬ 
main  without  further  approximations,  but  sometimes  a  finite  number  of  series 
expansion  terms  is  used  to  fit  the  linear  chromatic  dispersion  of  a  medium  or 
of  a  waveguide.  What  we  understand  under  NEE  in  the  following  assumes  an 
exact  treatment  of  the  dispersion  operator. 

Next,  we  approximate  the  nonlinear  coupling  term.  Unlike  in  NLS,  we 
preserve  the  frequency  dependence  exactly,  but  neglect  the  transverse  wave- 
number  dependence: 


Q(kx,ky,u)  — 


(ui  -  Wr)  +  WR 


2e0c2v/w2e(w)/c2-fc2-A;2  2  e0cn(wR) 


(1.56) 


Here,  as  in  the  free  propagation  term,  we  neglect  the  chromatic  dispersion  of 
the  background  index  of  refraction. 

After  putting  the  above  approximations  for  K  and  Q  into  the  original 
UPPE,  we  obtain 


dzA  =  ivgl{w  -  wr)A  +  iD(ui  -  ujr)A 


ic 


+ 


2wru(wr) 

WR 

2eocrc(u>R) 


(l  +  : 


WR 

(1  +  ^Z^)V 

WR 


'-)-\kl  +  kl)A 


(1.57) 


Finally,  transforming  into  the  real-space  representation,  we  arrive  at  NEE 

dzA  +  vg  ldtA  =  iD(idt)A+  xr- (1 4 — —dt)  1  A±A+  - — -fe — r(l  H — —dt)V 
9  2fcR  «r  2e0  ng(wR)  wR 

(1.58) 

Thus,  the  additional  approximations  underlying  the  NEE  are  paraxiality  both 
in  the  free  propagator  and  in  the  nonlinear  coupling,  and  a  small  error  in 
the  chromatic  dispersion  introduced  when  the  background  index  of  refraction 
is  replaced  by  a  constant,  frequency  independent  value  in  both  the  spatio- 
temporal  correction  term  and  in  the  nonlinear  coupling  term.  Note  that  the 
latter  approximations  are  usually  not  serious  at  all. 
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Fig.  1.1.  Supercontinuum  generated  in  a  femtosecond  pulse  propagating  in  air.  Pull 
curve  was  obtained  from  the  full  UPPE  simulation  while  the  other  curves  correspond 
to  NEE  equation  simulated  with  two  and  three  terms  included  in  the  dispersion 
operator. 


As  in  all  envelope  equations  a  reference  frequency  and  a  reference  wave- 
number  appear  in  the  NEE.  They  are  chosen  equal  to  the  central  frequency 
and  wave-number  of  the  input  pulse  in  practical  calculation,  but  one  has 
to  keep  in  mind  that  these  quantities  are  artificial  and  to  a  certain  de¬ 
gree  arbitrary  “gauge”  parameters  that,  of  course,  do  not  appear  in  the 
Maxwell’s  equations.  Consequently,  numerical  solutions  should  not  depend  on 
how  the  reference  is  chosen.  In  other  words  a  propagation  equations  should 
be  “reference-frequency-invariant.”  While  NLS  is  manifestly  dependent  on  the 
reference  choice,  NEE  is  nearly  invariant  although  wr  appears  in  it  several 
times.  For  example,  the  spatio-temporal  focusing  correction  term  (operator) 
WrJ(1  +  that  appears  in  the  NEE  equation  and  modifies  the  diffrac¬ 

tion  term  seems  to  depend  on  the  reference  wr,  but  it  is  in  fact  proportional 
to  w_1  as  long  as  it  is  implemented  in  the  spectral  domain  that  allows  to 
include  all  orders  of  the  series  expansion.  It  is  to  be  stressed  that  this  (ap¬ 
proximate)  invariance  is  only  achieved  in  the  infinite  order.  Truncating  the 
operators  that  appear  in  the  NEE  to  finite  number  of  series  expansion  terms 
breaks  the  invariance  and  brings  about  undesirable  artifacts  as  we  point  out 
in  the  following  example. 

We  consider  a  25-femtosecond  (0.1mm  waist)  pulse  with  a  carrier  wave¬ 
length  of  775  nm  and  power  of  8GW  propagating  in  air.  The  pulse  duration  is 
chosen  very  short  to  highlight  propagation  effects  that  are  absent  in  the  NLS 
approach,  namely  space-time  focusing  and  the  frequency  dependent  nonlinear 
response  (shock  formation).  We  compare  supercontinuum  generation  using 
UPPE  with  full  chromatic  dispersion  of  dry  air  taken  into  account  in  the 
wavelength  region  from  1200  to  200  nm  [37]  and  the  NEE  equation  with  the 
dispersion  operator  T>  =  Yhk=2  /3(,k^ / kl(idt)k  (Eqn.  (??))  is  expanded  up  to 
the  second  ( L  =  2)  and/or  third  (L  =  3)  order  with  j3l'k>  being  purely  real  in 
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this  case.  We  term  the  latter  approximations  ad2NEE  and  ad3NEE,  respec¬ 
tively  (standing  for  NEE  with  approximate  dispersion).  In  all  cases  we  assume 
an  instantaneous  optical  Kerr  effect  with  n.2  =  5  x  10_23m2/W,  and  plasma 
generation  by  multiphoton  ionization.  The  reader  is  referred  to  Ref.  [28]  for  a 
physical  description  of  the  model. 

Figure  1.1  shows  the  pulse  spectrum  after  the  self-focusing  collapse  is 
arrested  by  plasma  generation.  In  all  cases,  a  broad  high-frequency  component 
is  generated  on  the  steepened  trailing  edge  of  the  pulse  as  described  previously 
in  Ref.  [38].  However  the  details  of  the  spectra  are  rather  different.  Here,  the 
UPPE  solution  describes  the  correct  propagation  properties  of  all  wavelengths 
that  contribute  to  the  spectral  range  shown.  The  difference  between  the  UPPE 
and  ad2(3)NEE  solutions  can  be  traced  to  difference  in  the  susceptibility  they 
model.  It  happens  that  the  GVD  is  rather  small  around  800  nm  and  the 
approximated  susceptibility  rapidly  deviates  from  the  actual  susceptibility  at 
higher  frequencies.  Including  the  third-order  dispersion  substantially  improves 
the  agreement  with  the  UPPE  solution.  The  remaining  discrepancy  is  then 
restricted  to  the  high-frequency  range  in  which  the  supercontinuum  spectral 
intensity  falls-off.  This  demonstrates  that  in  the  NEE  the  dispersion  operator 
should  be  treated  exactly  in  the  spectral  domain  or  care  should  be  exercised 
in  approximating  chromatic  dispersion  by  an  expansion.  When  the  dispersion 
is  handled  properly,  NEE  is  an  excellent  approximation.  It  can  be  shown  that 
the  error  it  introduces  is  of  fourth  order  in  the  transverse  wave-number. 

1.3.3  Partially  corrected  NLS 

The  Partially  Corrected  NLS  (PC-NLS)  equation  can  be  viewed  as  a  “sim¬ 
plification”  of  NEE.  It  is  derived  from  the  UPPE  in  the  same  way,  with  one 
additional  step.  Namely,  the  following  first  order  series  expansion  is  applied 
in  the  correction  term  of  the  free  propagator  in  Eqn.(1.57): 

+  (1.59) 

wr  wr 

This  step  is  meant  to  make  it  easy  to  implement  a  numerical  solver  in  the  real 
space,  as  it  results  in  the  equation  that  only  contains  “simple”  differential 
operators  in  the  real-space  representation: 

dzA  +  v~ldtA  =  iD(idt)A  +  -i-(l  -  —  dt)A±A  +  ■?*  -.-(I  +  —dt)V 

9  2kRs  uR  2  e0n£(wR)  wR 

(1.60) 

While  it  may  seem  that  the  Partially  Corrected  NLS  is  essentially  NEE 
with  a  “little  more”  approximation,  this  equation  is  not  to  be  recommended. 
Because  of  the  arbitrary  truncation  of  an  infinite  series,  the  dispersion  prop¬ 
erties  of  the  linear  part  of  this  equation  are  unphysical.  While  the  PC-NLS 
provides  better-than-NLS  approximation  around  the  reference  frequency  wR, 
its  dispersion  properties  become  rather  pathological  around  w  ss  2 wR  where 
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Fig.  1.2.  Spatial  spectra  of  a  supercontinuum  generated  in  air.  Left  panel  repre¬ 
sents  a  solution  obtained  from  UPPE,  the  right  panels  is  a  corresponding  spectrum 
obtained  from  PCNLS.  Artifacts  in  the  spectrum  around  the  wavenumbers  that 
correspond  to  the  twice  the  reference  frequency  are  clearly  visible. 


its  diffraction  term  changes  sign  as  a  consequence  of  the  truncated  correction 
factor.  Artifacts  in  the  angular  distribution  of  the  spectrum  can  be  observed  at 
high  frequencies  beyond  u>  «  2(Vr.  This  is  illustrated  in  Fig.  1.2.  Consequently, 
this  equation  is  only  applicable  in  the  same  regime  as  the  NLS,  namely  when 
the  spectrum  of  the  pulse  remains  relatively  narrow. 


1.4  Applications  of  UPPE  simulators 

This  section  provides  three  illustrative  applications  of  the  z-UPPE  model. 
The  first  is  the  computationally  more  challenging  as  it  involves  a  full  3D  + 
time  simulation  of  the  propagation  of  a  wide  pancake  shaped  pulse  in  air. The 
second  provides  a  nice  illustration  of  the  need  to  go  beyond  the  paraxial 
approximation  for  nonlinear  X-wave  generation  in  condensed  media  and  the 
last  illustrates  the  subtle  interplay  between  plasma  generation  and  chromatic 
dispersion  in  limiting  the  extent  of  the  supercontinuum  spectrum. 

1.4.1  Femtosecond  Atmospheric  Light  Strings 
Multiple  filament  formation 

This  application  of  the  z-UPPE  simulator  illustrates  break-up  of  a  high-power, 
wide  femtosecond  pulse  into  chaotically  interacting  light  filaments.  Beams  that 
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carry  power  exceeding  the  critical  self-focusing  power  many  times  usually 
break  transversally  into  multiple  filaments.  To  capture  such  dynamics,  a  fully 
spatially  resolved  simulator  is  needed  that  doesn’t  impose  axial  symmetry. 

Below,  we  illustrate  how  such  multiple  filaments  are  concurrently  created 
at  different  transverse  and  longitudinal  locations,  and  how  they  interact  with 
the  low-intensity  background.  It  was  first  proposed  in  [15]  that  such  an  in¬ 
teraction  is  crucial  for  long  distance  propagation  of  high-power  femtosecond 
pulses  in  air.  The  basic  idea  is  that  of  dynamic  exchange  of  energy  between 
multiple,  essentially  unsynchronized  and  spatially  sharply  localized  filament 
cores  and  the  low-intensity,  spatially  wide  pedestal  of  the  beam. 

In  this  wide-beam  simulation,  the  initial  condition  is  a  Gaussian  pulse 
with  a  phase  perturbation.  The  waist  of  the  initially  collimated  Gaussian  was 
chosen  to  be  5mm,  the  pulse  duration  is  500  fs,  A  =  248nm,  and  the  maximal 
intensity  is  2  x  1014Wm~2.  The  total  pulse  energy  is  approximately  9mJ.  A 
random  phase  perturbation  is  imposed  on  the  pulse  to  initiate  the  transverse 
break-up  of  the  pulse  into  multiple  filaments  (see  1.3:).  We  adjusted  the 
amplitude  of  the  perturbation  such  that  it  results  in  the  filamentation  onset 
after  a  few  meters  of  propagation. 

The  first  stage  of  nonlinear  self-focusing  is  driven  by  the  smooth,  large- 
scale  profile  of  the  pulse.  After  a  few  meters,  local  perturbations  develop  into 
hot  spots  which  grow  into  high-intensity  filaments.  The  first  panel  shows  the 
overall  scale  of  the  input  pulse  with  the  high-intensity  regions  forming  from 
the  low- intensity  background.  There  is  practically  no  plasma  formation  at  this 
propagation  distance. 

The  initial  perturbations  grow  rapidly  and  reach  intensities  high  enough 
to  ionize  air  (second  panel).  Collapse  of  a  filament  is  eventually  regularized 
by  plasma  induced  defocusing.  That  causes  decay  of  the  filament  and  returns 
most  of  its  energy  into  the  low-intensity  background.  Prom  there,  new  fil¬ 
aments  grow  and  these  replenishment  cycles  repeat  with  relatively  modest 
energy  losses  to  plasma  generation  (subsequent  panels). 

Later  in  the  propagation,  filaments  start  to  appear  in  the  peripheral  re¬ 
gions  further  from  the  center.  This  is  due  to  less  overall  intensity  and  therefore 
slower  self-focusing  and  growth  of  perturbations.  Though  it  is  not  evident  on 
these  fluence  pictures,  later-stage  filaments  tend  to  generate  less  plasma  than 
the  ones  that  appear  at  the  very  beginning  of  the  filamentation  onset.  This 
is  the  stage  when  the  single-filament  dynamic  spatial  replenishment  scenario 
crosses  over  to  a  regime  where  replenishment  energy  originates  in  “neighbor¬ 
ing”  filaments  rather  than  from  the  same  one. 

This  gradual  change  in  the  dynamics  is  reflected  also  in  the  plasma  gener¬ 
ation  as  shown  in  Fig.  1.4  Initial  sharp  spikes  in  the  total  number  of  generated 
electrons,  associated  with  the  onset  of  individual  filaments,  decay  with  dis¬ 
tance  We  expect  the  shot-to-shot  fluctuation  to  smooth-out  the  sharp  features 
of  this  curve  due  to  randomization  of  the  filament  formation.  The  late-stage 
filaments  are  less  “organized”  than  those  created  just  after  the  self-focusing 
onset.  Consequently,  it  takes  less  of  the  plasma  generation  to  arrest  their  col- 
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Fig.  1.3.  Transverse  profile  of  the  fluence  (time-integrated  intensity)  “measured” 
at  several  propagation  distances.  The  color  scheme  was  chosen  such  that  one  can 
see  the  structure  of  the  low-intensity  background.  On  the  other  hand,  it  makes  it 
difficult  to  compare  filament  intensities.  The  size  of  the  depicted  domain  in  all  panels 
is  1cm  x  lcm. 


lapse.  One  can  say  that  the  increasing  “disorder”  in  the  developing  composite 
pulse  makes  the  collapse  arrest  due  to  plasma  more  efficient  and  thus  con¬ 
tributes  to  the  ability  of  the  pulse  to  propagate  over  long  distances.  One  can 
speculate,  and  recent  experiments  indicate  that  a  regime  can  be  eventually 
reached  where  the  plasma  generation  is  almost  negligible. 
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Fig.  1.4.  Total  number  of  free  electrons  generated  as  a  function  of  propagation 
distance. 


1.4.2  Nonlinear  X-waves  in  Condensed  Media 

In  this  Section,  we  present  an  example  of  simulation  of  an  ultrashort  pulse  in 
water.  This  particular  example  requires  that  the  solver  captures  correctly  the 
linear  dispersion  in  a  broad  range  of  frequencies  and  propagation  angles,  and 
is  thus  an  ideal  candidate  for  UPPE  application.  Full  details  can  be  found  in 
an  earlier  publication  [39] 

We  consider  a  loosely  focused  femtosecond  pulse  centered  around  the 
520  nm  wavelength,  propagating  in  a  water  sample.  An  appropriate  combina¬ 
tion  of  focusing  and  pulse  intensity  and  duration  results  in  a  long  (compared 
to  the  Rayleigh  range  corresponding  to  the  transverse  size  of  the  beam  at  the 
water-cell  entrance) . 

Figure  1.5  shows  the  resulting  filament  transverse  dimension  (size)  for 
several  pulse  energies  as  a  function  of  the  propagation  distance  in  water  (left). 
The  right  panel  illustrates  that  supercontinuum  is  easily  generated  at  these 
energies.  We  thus  deal  with  a  highly  non-NLS  regime. 

The  question  we  want  to  shed  light  on  in  this  numerical  experiment  is  what 
mechanism  is  responsible  for  creation  of  that  seemingly  several  centimeters 
long  filament.  Further,  we  want  to  know  if  the  mechanism  is  universal  in  any 
way. 

First,  it  is  important  to  note  that  what  is  actually  observed,  in  exper¬ 
iment  and  in  simulation  alike,  is  not  a  “steady-state”  self-guided  filament. 
Rather,  we  deal  with  a  series  of  pulse  splitting  events  akin  to  the  scenario  of 
spatial  dynamical  replenishment  [7].  In  this  case,  however,  the  role  of  plasma 
as  the  arrestor  of  the  self-focusing  collapse  is  less  pronounced  compared  to 
propagation  in  air. 
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Fig.  1.5.  Femtosecond  pulse  loosely  focused  into  water  creates  a  long  filament 
with  a  nearly  constant  diameter  that  extends  over  many  Rayleigh  ranges  (left), 
and  generates  a  broad  supercontinuum  spectrum  (right).  This  seemingly  stationary 
filament  is  created  by  a  series  of  very  dynamic  multiple  pulse  splittings  illustrated 
in  the  following. 


Fig.  1.6.  Pulse  splitting  cycle:  Single  central  pulse  undergoes  splitting.  The  split- 
off-pulses  act  a  “scatterers,”  which  concentrate  most  of  the  energy  in  the  spatio- 
temporal  spectrum  around  loci  that  support  “diffractionless”  wave-forms.  This  pro¬ 
vides  the  energy  for  a  new  central  pulse,  and  the  cycle  repeats. . . 


Figure  1.6  shows  a  series  of  snapshots  that  depict  the  temporal  profile  of 
on-axis  intensity  of  the  now  quite  complicated  “pulse.”  One  observes  several 
cycles  consisting  of  formation  and  subseqent  splitting  of  a  sub-pulse  in  the 
center  of  the  time  domain.  The  “daughter”  sub-pulses  resulting  from  each 
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pulse-split  event  play  an  important  role  in  the  formation  of  the  spatial  and 
temporal  spectrum.  Namely  they  are  still  intense  enough  to  induce  localized 
changes  to  the  material  susceptibility  that  in  turn  follow  these  split-off  pulses 
and  thus  propagate  with  different  “group”  velocities.  These  material  waves 
then  act  as  scatterers  in  a  three-wave  mixing  process  that  transforms  the  input 
optical  waves  into  scattered  ones.  Linear  propagation  dispersion  properties 
together  with  the  propagation  velocity  of  the  material  wave  then  determine 
where  in  the  spectral  space  will  the  scattered  energy  accummulate. 
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Fig.  1.7.  Angle-resolved  spectrum  of  an  ultrashort  pulse  propagating  in  water  (left 
panel).  The  dashed  lines  represent  the  loci  of  spectral  energy  concentration  predicted 
from  an  effective  three-wave  mixing  argument.  The  right  panel  compares  these  loci 
to  the  manifold  that  supports  the  spectrum  of  z-invariant  X-wave  solutions  that 
propagate  without  distortions. 


The  resulting  spatio-temporal  spectrum  of  a  loosely  focused  ultrashort 
pulse  after  propagation  in  water  is  shown  in  Fig.  1.7  on  left.  The  dashed  lines 
(in  both,  left  and  right  panels)  represent  the  loci  where  energy  concentrates 
due  to  the  non-linear  interactions  irrespectively  of  the  details  of  the  underly¬ 
ing  dynamics.  The  resulting  central  X-shaped  feature  is  always  close  to  the 
manifold  (shown  in  full  line,  right  panel)  that  supports  the  z-invariant  X- 
waves  that  propagate  long  distances  without  changing  their  spatio-temporal 
shapes.  In  any  normal-GVD  medium,  the  “theoretical  X-wave”  spectrum  and 
the  “real-pulse”  spectral  concentration  will  be  close  to  each  other  because  of 
the  simple  landscape  of  chromatic  dispersion  in  the  space  of  frequency  and 
transverse  wavenumber.  Thus,  even  highly  non- stationary  pulses  inherit  their 
tendency  for  long-distance  propagation  form  the  nonlinear  X-waves. 

1.4.3  Supercontinuum  Genration  in  Bulk  Media 

In  this  Section,  we  describe  comparative  simulations  that  provided  further 
insight  into  how  supercontinuum  is  generated  by  powerful  femtosecond  pulses 
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propagating  in  bulk  medium.  Readers  interested  in  more  details  can  found 
those  in  [22]  and  in  [21] 

The  SC  generation  requires  high  intensities  that  are  initiated  by  self- 
focusing  collapse  in  the  medium.  The  mechanism  that  arrests  the  collapse 
is  supplied  by  multi-photon  ionization  (MPI)  which  produces  both  a  di¬ 
rect  energy  loss  for  the  collapsing  field  and  an  electron-ion  plasma,  which 
subsequently  absorbs,  defocuses,  and  spectrally  blue-shifts  the  optical  field. 
The  combined  effects  of  MPI  and  plasma  defocusing  provide  a  mechanism 
for  the  arrest  of  self-focusing  collapse,  which  clamps  the  maximum  inten¬ 
sity  Imax  reached  by  the  collapsing  pulse,  and  also  limits  the  maximum 
rate  of  plasma  generation  dtp  oc  with  p  the  plasma  density  and  K 

the  order  of  the  MPI.  The  standard  scenario  [23,40-43]  maintains  that  the 
dominant  contribution  to  spectral  broadening  comes  from  the  spectral  blue- 
shifting  due  to  the  plasma  [42,44],  with  the  maximum  blue-side  frequency 
shift  Aumax  oc  dtp  oc  ImaX-  Thus,  according  to  the  standard  scenario,  limit¬ 
ing  of  the  maximum  intensity  Imax  as  self-focusing  is  arrested  is  the  dominant 
factor  that  determines  the  spectral  extent  Au)max  of  SC  generation. 


Fig.  1.8.  “Artificial”  vs  “real”  water  susceptibility  used  in  comparative  simulations 
designed  to  test  the  standard  supercontinuum  generation  scenario  for  bulk  media. 


In  order  to  put  the  standard  SC  generation  scenario  to  a  test,  we  compare 
the  water  SC  generation  simulations  with  analogous  simulations  performed 
using  an  artificial  medium  which  is  the  same  as  the  “original”  but  with  a 
modified  linear  dispersion.  The  later  is  “constructed”  such  that  the  artificial 
medium  exhibits  self-focusing  and  plasma  dynamics  that  are  almost  identical 
to  those  of  the  real  medium  model. 

Frequency  mapping  ui  — ►  J?(w)  is  used  to  construct  the  artificial  water 
susceptibility  is  shown  in  long-dashed  line  in  Fig.  1.8  (the  dotted  line  shows 
identity  function  for  comparison ).  The  artificial  and  water  susceptibility  func¬ 
tions  are  depicted  in  dashed  and  full  lines,  respectively  in  Fig.  1.8.  The  arrow 


24  J.  V.  Moloney,  M.  Kolesik:  Ultrashort  Pulse  Propagators 

illustrates  the  susceptibility  transformation.  Susceptibility  is  preserved  around 
the  pulse  carrier  frequency,  in  order  to  preserve  all  the  quantities  that  control 
the  SC  generation  under  the  standard  scenario,  such  as  the  highest  light  inten¬ 
sities  and  plasma  density.  Consequently,  according  to  the  standard  scenario, 
the  two  media  should  generate  the  same  SC  spectra.  Numerical  simulations 
are  used  to  check  if  this  is  the  case. 


Fig.  1.9.  Local  susceptibility  variation  induced  by  the  nonlinear  effects  in  the  fem¬ 
tosecond  pulse.  Water  and  a  the  comparative  artificial  medium  exhibit  very  similar 
responses. 


To  show  that  both  model  media  should  indeed  provide  essentially  the  same 
SC  spectrum  in  our  simulations,  provided  the  standard  SC  scenario  is  correct, 
we  present  a  comparison  of  the  nonlinear  response  for  water  and  the  artificial 
medium  in  Fig.  1.9.  The  nonlinear  response  consist  in  a  local  change  A\  of 
the  medium  susceptibility.  As  a  function  of  the  “local  time”  at  pulse  location, 
it  first  increases  due  to  the  increasing  intensity  created  by  the  selffocusing. 
Then  the  susceptibility  change  decrease  into  negative  values  because  of  the  de- 
focusing  caused  by  the  plasma  generated  in  the  high-intensity  light  pulse.  The 
picture  shows,  that  the  nonlinear  response  of  the  artificial  medium  is  extremely 
close  to  that  of  the  water.  This  indicates  that  if  the  standard  model  of  SC 
generation  is  complete,  the  almost  same  response  dynamics  implies  closely 
similar  SC  spectra. 

The  comparison  of  the  SC  spectra  obtained  using  the  realistic  and  arti¬ 
ficial  water  susceptibilities  is  shown  in  Fig.  1.10.  The  original  and  artificial 
medium  spectra  agree  quite  well  in  the  vicinity  of  the  excitation  wavelength. 
However,  at  high  frequencies  the  two  materials  produce  drastically  different 
continua.  The  dashed  line  shows  the  SC  spectrum  transformed  by  the  same 
transformation  that  produced  the  artificial  medium  susceptibility  from  the 
water  susceptibility.  This  appears  to  be  very  close  to  the  artificial  medium 
spectrum.  This  indicates  that  the  extend  of  the  spectrum  is  actually  deter¬ 
mined  mostly  by  the  linear  dispersion  properties  of  the  medium. 
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Fig.  1.10.  Supercontinuum  spectra  generated  in  water  and  in  the  artificial  medium. 
Dashed  line  represents  the  spectrum  edge  obtained  from  the  real  water  spectrum 
the  same  way  as  the  artificial  susceptibility  was  obtained  from  the  original  water 
susceptibility.  The  standard  SC  scenario  predicts  the  same  spectra  for  both  models. 


Thus,  we  arrive  at  the  conclusion  that  the  standard  scenario  for  supercon- 
tinuum  generation  in  bulk  media  is  incomplete.  Although  it  correctly  identifies 
the  key  quantities  (  peak  intensity,  plasma  density  generation  rate)  and  pro¬ 
cesses  (collapse  arrest,  MPI),  our  numerical  experiments  demonstrate  that  it 
doesn’t  explain  the  supercontinuum  spectral  properties. 

To  complete  the  supercontinuum  generation  picture  for  bulk  media,  the 
medium’s  linear  susceptibility,  as  a  function  of  frequency,  must  be  taken  into 
account  as  a  key  factor  that  determines  the  maximal  attainable  width  of  the 
white  light  spectrum. 

Further  simulations  (not  shown)  also  showed  that  the  so-called  supercon¬ 
tinuum  band-gap  dependence  is  due  to  stronger  chromatic  dispersion  at  higher 
frequencies.  Namely,  increasing  the  excitation  frequency  has  a  similar  effect 
on  the  SC  spectrum  as  the  artificial  modification  of  the  susceptibility  we  used 
in  our  simulations;  The  curvature  of  the  dispersion  function  is  stronger  closer 
to  the  band-gap  and  reduces  the  spectral  broadening  due  to  the  larger  phase- 
mismatch  in  the  underlying  wave  mixing  processes. 


1.5  Future  Work  and  Perspectives 

This  chapter  has  focused  on  procedures  for  deriving  ultra  short  propagation 
pulse  propagation  models  starting  from  Maxwell’s  equations.  The  unidirec¬ 
tional  pulse  propagation  equation(UPPE)  possesses  all  of  the  essential  ingre¬ 
dients  for  current  and  future  studies  in  the  emerging  field  of  extreme  non¬ 
linear  optics.  We  have  shown  that  the  latter  provides  a  unifying  theoretical 
framework  from  which  the  many  propagation  equations  in  the  literature  can 
be  seamlessly  derived.  Moreover,  the  physical  approximations  made  in  these 
models  are  clearly  exposed  and  their  interrelationships  stressed. 
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